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A fundamental mystery in earthquake physics is “how can an earthquake be triggered by distant 
seismic sources?” Here, we use discrete element method simulations of a granular layer, during 
stick-slip, that is subject to transient vibrational excitation to gain further insight into the physics 
of dynamic earthquake triggering. Using Coulomb friction law for grains interaction, we observe 
delayed triggering of slip in the granular gouge. We find that at a critical vibrational amplitude 
(strain) there is an abrupt transition from negligible time-advanced slip (clock advance) to full 
clock advance, i.e ., transient vibration and triggered slip are simultaneous. The critical strain is 
order of 10 —6 , similar to observations in the laboratory and in Earth. The transition is related to 
frictional weakening of the granular layer due to a dramatic decrease in coordination number and the 
weakening of the contact force network. Associated with this frictional weakening is a pronounced 
decrease in the elastic modulus of the layer. The study has important implications for mechanisms 
of triggered earthquakes and induced seismic events and points out the underlying processes in 
response of the fault gouge to dynamic transient stresses. 
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I. INTRODUCTION 

Dynamic triggering of earthquakes by seismic waves is 
a robustly observed phenomenon that is well-documented 
for over 30 major earthquakes worldwide [f| and many 
more smaller earthquakes 00. Recent observations 
based on new, more sensitive instrumentation show that 
a majority of earthquakes may be dynamically trig¬ 
gered 0I3- Laboratory-scale experiments and seismo- 
logical observations indicate that a key role in dynamic 
earthquake triggering may be played by granular mate¬ 
rials, termed “fault gouge”, accumulated at the core of a 
geologic fault [8, 9J. The observations at the laboratory 
and field scales strongly suggest that the nonlinear dy¬ 
namical response of the gouge material significantly con¬ 
tributes to triggering, although details remain unquanti¬ 
fied. Direct access to the earthquake fault gouge with¬ 
out changing its microstructure and loading history is 
not possible. However, we here aim to characterize the 
granular physics of triggering on laboratory scales using 
physical experiments and numerical simulations. 

Granular layers exhibit stick-slip dynamics when they 
are subjected to shearing, at sufficiently high confin¬ 
ing pressures and low shearing velocities [IMS]. The 
stick-slip instabilities have been associated with a non¬ 
monotonic shear stress vs shear strain response of gran¬ 
ular materials that have frictional constituents or fric¬ 
tional dissipation m- The stick-slip dynamics are anal¬ 
ogous to the seismic cycle in earthquake fault systems. 
Fault systems accumulate strain energy during the in- 
terseismic period of the seismic cycle, just as a sheared 
granular layer does during the stick phase of the stick- 
slip cycle mun]. Laboratory scale observations confirm 
that mechanical vibrations with adequate amplitudes can 
change the mechanical and frictional properties of the 
granular layer, changing its macro-scale response. This 
includes a transition from a solid-like behavior to a tran¬ 
sient, fluid-like one [T8jf24j . The behavior of granular 
materials under different loading conditions and to dif¬ 
ferent perturbations is controlled by their evolving inter¬ 
nal structure including the contact force networks, par¬ 
ticle rearrangements and force distribution between the 
particles inside the granular layer f25ll27l . Jia et al |23] 
identified two regimes of fast nonlinear dynamics versus 
the input amplitude. In the first regime, the interac¬ 
tion between sound waves and the granular medium is 
reversible: during the wave excitation the modulus can 
change. However, neither velocity nor sample density are 
changed after the wave passage and the force network re¬ 
mains nearly unchanged. In the second regime, beyond 
a certain amplitude threshold depending on the applied 
load, the sound-matter interaction becomes irreversible. 
In addition, the wave velocity and corresponding elastic 
modulus remain weakened after the wave transient, and 
permanent deformation is observed corresponding to an 
accompanying compaction. This finding highlights the 
relationship between the macroscopic elastic weakening 
and the local change of the contact network, induced by 


strong sound vibration in the absence of visible grain mo¬ 
tion [23] . 

Johnson et al. [8] observed both instantaneous and 
delayed triggered (cascading) slip in the lab, when vi¬ 
bration amplitudes corresponding to strains 10 -6 
are applied at shear stress levels of « 95% of the fail¬ 
ure value. Other studies also demonstrated the exis¬ 
tence of threshold values of strain amplitude for dynamic 
earthquake triggering. The existence of a unique strain 
threshold value is an important open question [6], how¬ 
ever there is increasing evidence that in many cases dy¬ 
namic earthquake triggering may be governed by such 
a threshold mechanism 0 |U HH]. Johnson et al. 0 
also observed other features in common with earth faults 
including disruption in the earthquake recurrence inter¬ 
val (the time interval between earthquakes) in response 
to dynamic perturbations, as well as triggering-induced 
changes in the gouge material modulus. We have previ¬ 
ously investigated the deformation characteristics for dy¬ 
namically triggered slip and the influences of vibration 
amplitude on triggering using two dimensional discrete 
element simulations of a granular fault gouge [29, [30]. 
In a follow-up study, we developed a three dimensional 
discrete element model of the granular fault gouge that 
showed similar dynamics to earthquake stick-slip cycles 
and aseismic creep prior to slip events [31, =22] • In that 
system, we have characterized the short and long term 
influences of triggering with regard to stick-slip size dis¬ 
tribution and recurrence intervals[33l [34] . Here, we re¬ 
port results of three dimensional discrete element method 
simulations of granular gouge layers subjected to bound¬ 
ary vibrations to understand the grain-scale mechanics 
of dynamic earthquake triggering and the existence of 
a triggering threshold under certain confining pressures. 
We also explore the nature of the transition to significant 
clock-advanced triggered slip events at different trigger¬ 
ing vibration amplitudes and frequencies. This modeling 
work is set to complement experimental observations ob¬ 
tained from double-direct shear experiments by Johnson 
et al. [8]. 


II. MODEL SETUP 

Figure [l] illustrates the simulated granular gouge layer. 
The model consists of three layers of particles: a driv¬ 
ing block at the top, a granular gouge layer and a sub¬ 
strate block at the bottom. The driving and substrate 
blocks are used to confine the granular gouge by applying 
a constant normal force in the E-direction. The top driv¬ 
ing block moves at constant velocity in the positive X- 
direction and applies a shear force to the granular gouge 
layer. Each variable/parameter in our 3D DEM model 
is expressed in terms of the following basic dimensional 
units: L 0 = 150 /xm, to = 1 s and M 0 = 1 kg , for length, 
time and mass, respectively. Lq represents the largest 
particle radius within the overall DEM model. We run 
sheared granular layer simulations at a confining pres- 
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Numerical setup 


FIG. 1. 3D DEM model comprised of the driving block (top), 
a granular gouge layer (center) and a substrate block (bottom) 

sure of <r n = 6000 2 (40 MPa) and shearing velocity of 

Vx,o = 0.004^ ( 0.6 ^P) to achieve stick-slip dynamics. 
We chose parameter values to match laboratory experi¬ 
ments, rather than tectonic fault zones, although the two 
overlap in many ways. Further details about the model 
are provided in the supplementary materials. 

The simulations that are not subjected to vibration 
are called “reference” runs, while those with vibration 
are called “perturbed” runs. In the perturbed (also 
called “triggered”) simulation runs, an additional bound¬ 
ary condition consists of imposing a cyclic displacement 
in the Y direction for the bottom particles of the sub¬ 
strate for the duration of about At = 0.1 to- The char¬ 
acteristic of the vibration signal is shown in figure [ 2 }a 
(gray zone). The temporal displacement of the bound¬ 
ary vibration is described in the supplementary materi¬ 
als. Perturbations with longer duration have a bigger in¬ 
fluence than those with shorter duration. The influence 
of vibration duration is provided in the supplementary 
materials. Vibration normal to the boundary (displace¬ 
ment in Y direction) is found to be more effective in 
triggering slip than horizontal vibration (in X and Z di¬ 
rections) primalry due to better transmission of normal 
vibration to and through the cohesionless granular gouge 
layer compared with the horizontal vibration. The com¬ 
parison of horizontal and vertical vibration influences is 
included in the supplementary materials. For the DEM 
model, we use a range of vibration amplitudes includ¬ 
ing A ={ 1 ,10, 20,30,40, 50,60, 70,80,90,100} x Hr 7 L 0 . 
We can estimate the strain induced by the vibration as, 
e = where e is the induced strain, A is the vibration 

^vib 

amplitude in [Lq\ unit, and A^ is the vibration wave¬ 
length. The wavelength is A = -A-, where v is the sound 

J u % b 

speed and f v n, is the vibration frequency [35 . 

The reference 3D DEM model behavior are described 
in an earlier work m- The packing fraction in the sim¬ 
ulation is ^ 0.58 at the beginning of stick phases. The 
packing fraction gradually decreases while the granular 
layer dilates during the stick-phase to ~ 0.56. The stick- 
slip behavior is monitored by its friction coefficient time 
series. The friction coefficient, /i, is defined as the ratio 
of shear stress developed at the boundary layers to the 
imposed normal stress. 


III. RESULTS 

n our DEM model as well as the experimental setup 
studied by Johnson et al. [ 8 ], vibration amplitudes that 
induce a strain value of order ~ 10 -6 cause a time- 
advanced slip (clock advance) [33] . In both setups, the 
vibration causes an immediate weakening (reduction of 
shear strength) of the granular layer 

Figure [2}a shows the behavior for selected vibration 
amplitudes in the DEM model. Here, vibration ampli¬ 
tude larger than ~ 6 x 10 _6 Lo (corresponding to induced 
strain of ~ 3.9 x 10 -6 ) causes a sharp clock advance. 
Johnson et al. [ 8 ] report a highly perturbed stick-slip re¬ 
currence interval due to acoustic excitations, compared 
to the reference case. Similarly in the simulation, vibra¬ 
tion induced clock advance means that the recurrence 
time for the next event will be longer. In addition, John¬ 
son et al. [[36; report a similar vibration strain for induc¬ 
ing slow slip event in sheared granular fault gouge. We 
measure shear elastic modulus of the granular layer to 
monitor the evolution of elastic properties correspond¬ 
ing to the observed clock advance. The shear modulus 
of the granular layer is determined by applying a small 
shear strain cycle to the system both in the reference 
and in the perturbed simulations. For this purpose, the 
shearing of the granular layer is stopped and the system 
is given 10000 numerical time steps to relax from the 
shearing influences before performing the shear modulus 
measurements. We ensured that the average normal and 
tangential contact forces remain constant following the 
relaxation process. In addition, there have been no ob¬ 
servation of any failure event in the shear and normal 
stress signals and therefore no significant change in the 
contact network of the granular layer during this pro¬ 
cess. We then applied a cyclic shear strain to the top 
of the granular layer. We investigated the influence of 
different cyclic shear strain amplitudes for the modulus 
measurements. The results indicate that when the max¬ 
imum applied shear strain 7 maa , < 9 x 10 —6 , the shear 
modulus measurements are similar to each other. When 
the applied shear strain is further increased, the behav¬ 
ior becomes nonlinear. This observation suggests that 
shear strain values 7 max > 1 x 10 -5 induce consider¬ 
able particle contact rearrangement in the granular gouge 
layer. We therefore performed the shear modulus mea¬ 
surements with a load cycle of maximum strain ampli¬ 
tude 7 max ~ 8.6 x 10 -6 . Furthermore, the amplitude 
and duration of the protocol is selected such that it does 
not change contacts population during its application. 
A similar loading convention for modulus measurements 
has been suggested before by Nguyen et al. [37] . The 
shear modulus is determined by fitting a line to the initial 
unloading part of the stress-strain curve. The procedure 
is explained in the supplementary materials. 

Figure [2}b shows the shear modulus measurements 
for the reference and the perturbed simulations over the 
range of vibration amplitudes. During vibration, we ob¬ 
serve a decrease in the shear elastic modulus of the gran- 
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FIG. 2. (a) Friction coefficient signal, (b) shear elastic modulus measurements, (c) coordination number variations for reference 
(black) and perturbed simulations with a range of vibration strains. The vibration interval is shaded in all panels, subfigures 
(d), (e), and (f) show the variation of total number of contacts, number of weak contacts and number of strong contacts, 
respectively, during the full vibration interval. 


ular layer, indicating the mobilization of particles during 
the vibration interval. As the vibration amplitude in¬ 
creases, the shear modulus further decreases during the 
vibration interval. If the vibration amplitude is small 
(for this example, e < 3.9 x 10 -6 ), the shear modu¬ 
lus recovers after the vibration terminates. However, if 
e > 3.9 x 10 —6 , the shear modulus does not recover at 
the end of the vibration interval, but decreases further 
as the shearing continues beyond the vibration interval. 
This behavior leads to the clock advance of the large ex¬ 
pected slip event. This observation supports the hypoth¬ 
esis of Johnson and Jia [9] that, upon the application of 
sufficiently large vibration strains, the shear modulus of 
the granular medium decreases and the slip event occurs 
with a clock advance compared to its reference time. 

The variations of the coordination number ( CN ), the 
number of contacts per particle, for reference and se¬ 
lected perturbed runs are shown in Figure [2}c. The CN 
decreases slightly in the reference run as we approach 
the slip onset. During slip, the CN drops significantly. 
In the case of perturbed runs, the CN decreases signifi¬ 
cantly during the vibration interval. It recovers after vi¬ 
bration terminates to the reference level if the vibration 


strain is e < 3.9 x 10 -6 . However, for vibration strains 
e > 3.9 x 10 —6 , the CN initially recovers during the grad¬ 
ual removal of vibration, but it begins to decrease again 
leading to the clock advance of the slip event. Here, we 
note that the decreasing shear modulus of the granular 
gouge layer at the time of (triggered) slip as well as low 
values of CN are reminiscent of the response of an un¬ 
jammed granular layer. However, we cannot determine 
whether the system is jammed or shear-jammed [38] dur¬ 
ing stick phases because of the continuous slow rearrange¬ 
ment of the granular gouge layer that take place in the 
stick phase. Understanding the nature of the transition 
of the frictional granular layer from sticking to slipping 
will require further studies. 


IV. DISCUSSION 

The decrease of the CN despite the initial recovery is 
in agreement with the change in elastic shear modulus of 
the granular gouge layer. We now investigate the evolu¬ 
tion of the contact network of the granular gouge layer. 
The change in the total number of contacts during the 
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vibration interval is shown in figure [2^d and indicates a 
decrease of the total number of contacts as the vibration 
intensifies. The grain contacts show a slow recovery to¬ 
ward the end of the vibration interval for all vibrational 
strains. 

We now investigate the variation of the number of 
weak and strong contacts in the granular gouge layer. 
Weak contacts are those contacts that are carrying nor¬ 
mal forces smaller than the average normal contact force 
of the granular gouge layer. Strong contacts are carry¬ 
ing normal forces larger or equal to the average normal 
contact force of the granular gouge layer. The contact 
force distribution in granular systems is approximately 
an exponential distribution m- Strong and weak con¬ 
tacts form two subnetworks of the granular contact net¬ 
work with complementary mechanical properties [39, 00] . 
Strong contacts form the majority of the force chains or 
force bearing structure of the granular medium whereas 
weak contacts primarily support the strong contacts, dis¬ 
tribute the force over the whole granular contact network 
and help mobilizing the overall granular medium. The 
number of weak and strong contacts during the vibra¬ 
tion interval are shown in figure [2}e and [2}f, respectively. 
Figure shows that for vibration strains e > 3.9 x 10 -6 , 
the weak contact number decreases dramatically at peak 
of vibration. The number of strong contacts on the con¬ 
trary increases at peak vibration for e > 3.9 x 10 -6 as 
can be observed in figure [2]- f. The weak contacts regen¬ 
erate and recover as the vibration is removed, while the 
number of strong contacts decreases for vibration strains 
e > 3.9 x 10 -6 . These observations suggest that the 
strong contact network is strengthened when the vibra¬ 
tion strain increases initially. This is more or less ex¬ 
pected, because we are applying an excess confining pres¬ 
sure due to the vibration. Indeed, we observe some small 
scale fluctuations in the strong contact number synchro¬ 
nized with the applied vibration. Following the peak in 
the vibration strain, we observe a weakening of the strong 
contact network. After peak vibration, the weak contact 
recovers and strengthens above the reference level. This 
fact could explain the occurrence of delayed dynamic trig¬ 
gering: the strong contact network population is perma¬ 
nently depleted of members by the vibration application. 
It implies that the system is more susceptible to grain- 
scale rearrangements in correspondence of the continu¬ 
ously applied shear load, since the granular medium is 
in a critical state close to failure. This result indicates 
that the large vibration strain weakens the strong con¬ 
tact network of the granular layer and by this causes a 
clock advance of the expected slip event. 

We further explore the nature of the transition to sig¬ 
nificant clock advance for a range of vibration strains 
and frequencies. Figure [3} a shows the normalized clock 
advance time, defined as the difference between the ref¬ 
erence and triggered slip times divided by the recurrence 
interval of the reference slip event, for a range of vibra¬ 
tion strains applied in ^ 70% of the stick-phase of three 
different stick-slip cycles. This three stick-slip cycles be- 




FIG. 3. (a) Normalized clock advance time of triggered slip, 
defined as the difference between the reference and triggered 
slip times divided by the recurrence interval of the reference 
slip event, for three different reference slip events shown with 
symbols (circle), (triangle) and (square) at a range of vibra¬ 
tion strains, (b) Normalized clock advance time of triggered 
slip for a range of vibration strains at different triggering fre¬ 
quencies. Inset: minimum vibration strain required for caus¬ 
ing a significant clock advance at different vibration frequen¬ 
cies and ratios of vibration wavelength to mean particle di¬ 
ameter in the gouge layer corresponding to these frequencies 
(top axis of the inset). 


long to the same numerical setup, but is chosen from 
different temporal points in the stick-slip time-series of 
the system. The clock advance vs. vibration strain plot 
shows characteristics of a first order transition that takes 
place at vibration strain e ^ 10 _6 for all studied events, 
however the amount of clock advance can be different 
for different slip events. Due to the discontinuous nature 
of the granular materials, the fault gouge layer is sen¬ 
sitive to the vibration frequency, i.e. the frequency at 
which a wave packet explores the medium, introducing a 
competition between rearrangement length scale and vi¬ 
bration wavelength. We know from seismotectonics that 
an earthquake fault rupture generates seismic waves with 
a spectrum of frequencies, which influences the granular 
gouge of nearby or distant fault zones. It is thus of rel¬ 
evance to understand the consequences of triggering at 
different vibration frequencies. This is shown in Fig. [3} 
b where we perturb the granular gouge layer at a range 
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of frequencies from 250 to 2000 Hz. The transition is a 
first order type process for different vibration frequencies. 
The vibration strain at which the transition occurs de¬ 
creases with increasing the vibration frequency to about 
f v ib = 1500 Hz. The inset of Fig. [3]-b shows the minimum 
vibration strain required for a significant clock advance 
versus the frequency of the applied vibration. This fig¬ 
ure shows that by increasing the frequency of vibration, 
the minimum vibration strain required for observing a 
clock advance of the expected large slip event gradually 
decreases until we reach a vibration frequency of ^800 
Hz. From this point on by further increasing the vibra¬ 
tion frequency the minimum amplitude for clock advance 
only decreases slightly and stays at the vibration ampli¬ 
tude of e > 4 x 10 -6 for the vibration frequencies of 1500 
and 2000 Hz. 

This trend can be explained by calculating the vibra¬ 
tion wavelength. For example, for the vibration fre¬ 
quency f vib = 1000 Hz, we have X vib = j~ h = = 

2 x 10 _4 m, where A v i b and v are the vibration wavelength 
and speed of sound in the granular layer, respectively. 
The grain diameter in the granular gouge layer is in the 
range of [1.05; 1.5] x 10 -4 m. Thus we find that the vibra¬ 
tion wavelength for the vibration frequency f v i b = 1000 
Hz is very close to the average size of the granular gouge 
particles. Further increase of the vibration frequency de¬ 
creases the vibration wavelength. The wavelength at fre¬ 
quencies in the range of 800 to 1000 Hz can thus explore 
the granular media at the grain scale and perturb the con¬ 
tact network of the medium more effectively. However, 
decreasing the vibration frequency increases the vibra¬ 
tion wavelength. As a result, the perturbing vibration 
at frequencies lower than 500 Hz, are at length scales 
larger than the grain size and cannot significantly per¬ 
turb the grain arrangements and the contact network of 
the granular layer. A similar behavior for the trigger¬ 
ing frequency has been observed by Savage and Marone 
(3X1 in laboratory experiments with sheared and dynam¬ 
ically perturbed granular layers, where they notice that 
upon increasing the triggering force frequency, the fric¬ 
tion coefficient at dynamic (perturbed) failure starts to 
decrease compared to the friction coefficient at reference 
(unperturbed) failure until a critical frequency is reached. 
By further increase of the frequency beyond that critical 
value, the ratio of the friction at dynamic failure to the 
friction at the reference failure saturates [4Tj . The vi¬ 
bration frequency f v i b = 1000 Hz used in this study is 
in a range that can effectively perturb the granular layer 
at the grain scale and change its contact network signif¬ 
icantly with short time exposures. 


V. SUMMARY 

We have investigated the influence of boundary vibra¬ 
tion on stick-slip dynamics of sheared granular layers us¬ 
ing discrete element simulations. The numerical work 
is meant to complement the experimental observation 
by Johnson et al. [8]. In both discrete element simu¬ 
lations and experiments, a range of vibration strains is 
used for slip triggering and it is shown that above a criti¬ 
cal amplitude of about ~ 10 —6 , vibration causes a signif¬ 
icant clock advance of large slip events. Numerical sim¬ 
ulation shows that vibration initially influences the net¬ 
work of weak contacts in the granular gouge layer. The 
weak contacts can be regenerated after removal of vibra¬ 
tion, whereas strong contacts remain weakened compared 
to their pre-vibration state. We link the observed clock 
advance for e 10“ 6 to a major decrease of the co¬ 
ordination number as well as weakening of the strong 
contact network of the granular gouge layer. We further 
found that clock advance of the triggered slip event is 
a first order phase transition process in correspondence 
with increasing vibration strain amplitude. We would 
like to emphasize that in the Earth, the exposure time 
to dynamic stresses, sound velocity in the fault gouge 
and fluid pressure are among the factors that determine 
the effectiveness of triggering. The triggering threshold, 
if exists in Earth and field, may depend on the type of 
interaction between grains, for example on the existence 
of cohesion, humidity and fluid pressure (Scuderi et al. 
[42] 5. as well as on the physical and chemical properties of 
contact asperities and gouge roughness. These are among 
the areas that need further research and investigation. 
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